1 Mo_E18_5_MesCells Samples

We did scRNAseq (10X version 3) of mouse endochondral bones cell populations (in a balanced/even representation way, by FACS isolation) in perinatal stages E18 (just before birth) and postnatal day 1 (PN1).

1.1 About input files

The files used as input for this analysis are the files coming out of the the cellranger count pipeline.

In this case I used cellranger v6.0.0.

library(Matrix)
library(dplyr)
library(Seurat)
library(patchwork)
library(ggplot2)
library(DoubletFinder)
library(BiocStyle)
library(dittoSeq)
library(RColorBrewer)
ggplotColours <- function(n = 6, h = c(0, 360) + 15){
  if ((diff(h) %% 360) < 1) h[2] <- h[2] - 360/n
  hcl(h = (seq(h[1], h[2], length = n)), c = 100, l = 65)
}

1.2 Loading data from 10X

Seurat is able to analyze data from 10X experiments processed using CellRanger v3.

I first read all the data and then create a seurat object with the gene expression

# Load the Mo_E18_5_MesCells dataset
cluster_dir <- "~/scratch/stromal_cells/1-cellranger_mapping/Mo_E18_5_MesCells/Mo_E18_5_MesCells"
Mo_E18_5_MesCells.data <- Read10X(data.dir = paste0(cluster_dir,"/outs/filtered_feature_bc_matrix/"))
# Initialize the Seurat object with the raw (non-normalized data).
Mo_E18_5_MesCells <- CreateSeuratObject(counts = Mo_E18_5_MesCells.data, 
                                        project = "Mo_E18_5_MesCells",
                                        min.cells = 3, min.features = 200)
rm(Mo_E18_5_MesCells.data)
saveRDS(Mo_E18_5_MesCells, "Mo_E18_5_MesCells_init.rds")
Mo_E18_5_MesCells <- readRDS("Mo_E18_5_MesCells_init.rds")

1.3 Standard pre-processing workflow

The steps below encompass the standard pre-processing workflow for scRNA-seq data in Seurat. These represent the selection and filtration of cells based on QC metrics, data normalization and scaling, and the detection of highly variable features.

1.3.1 QC and selecting cells for further analysis

Seurat allows you to easily explore QC metrics and filter cells based on any user-defined criteria. A few QC metrics commonly used by the community include

The number of unique genes detected in each cell. Low-quality cells or empty droplets will often have very few genes Cell doublets or multiplets may exhibit an aberrantly high gene count Similarly, the total number of molecules detected within a cell (correlates strongly with unique genes) The percentage of reads that map to the mitochondrial genome Low-quality / dying cells often exhibit extensive mitochondrial contamination We calculate mitochondrial QC metrics with the PercentageFeatureSet function, which calculates the percentage of counts originating from a set of features We use the set of all genes starting with MT- as a set of mitochondrial genes

# The [[ operator can add columns to object metadata. This is a great place to stash QC stats
Mo_E18_5_MesCells[["percent.mt"]] <- PercentageFeatureSet(Mo_E18_5_MesCells, pattern = "^mt-")
# Show QC metrics for the first 5 cells
head(Mo_E18_5_MesCells@meta.data, 10)
##                           orig.ident nCount_RNA nFeature_RNA percent.mt
## AAACCCAAGGGTAATT-1 Mo_E18_5_MesCells      16021         3854   1.547968
## AAACCCACAGTTAAAG-1 Mo_E18_5_MesCells       6176         2005   2.347798
## AAACCCACATCCGATA-1 Mo_E18_5_MesCells      20024         4683   1.588094
## AAACCCAGTATCGATC-1 Mo_E18_5_MesCells       4619         1552   2.056722
## AAACCCAGTCAAATCC-1 Mo_E18_5_MesCells       8706         2499   1.435791
## AAACCCAGTCCATACA-1 Mo_E18_5_MesCells      17951         4668   2.016601
## AAACCCAGTCCCTCAT-1 Mo_E18_5_MesCells      19091         4038   1.498088
## AAACCCATCCATTTGT-1 Mo_E18_5_MesCells      18644         4030   2.558464
## AAACCCATCTCGAGTA-1 Mo_E18_5_MesCells       5175         1751   1.700483
## AAACGAACACCCTATC-1 Mo_E18_5_MesCells      20061         4628   1.031853

In the example below, we visualize QC metrics, and use these to filter cells.

We filter cells that have unique feature counts over 6,000 or less than 400 We filter cells that have >10% mitochondrial counts

# Visualize QC metrics as a violin plot
VlnPlot(Mo_E18_5_MesCells, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
        group.by = "orig.ident", ncol = 3, pt.size = 0)

 plot(x = Mo_E18_5_MesCells@meta.data$nCount_RNA,log = "y",
            xlab = "RNA counts",ylab = "Features",
      y = Mo_E18_5_MesCells@meta.data$nFeature_RNA, 
      type = "p",pch=20 , cex= 0.2, col="darkblue")
abline(h = c(400,6000), lty=2)
abline(v = c(1000), lty=2)

 plot(x = Mo_E18_5_MesCells@meta.data$nFeature_RNA,
      xlab = "Features",ylab = "perc.mit",
      y = Mo_E18_5_MesCells@meta.data$percent.mt, 
      type = "p",pch=20 , cex= 0.2,#xlim = c(0,35000),
      col="darkblue")
abline(v = c(400,6000), lty=2)
abline(h = c(10), lty=2)

Mo_E18_5_MesCells <- subset(Mo_E18_5_MesCells, subset = nFeature_RNA > 400 & nFeature_RNA < 6000 & percent.mt < 10)
dim(Mo_E18_5_MesCells)
## [1] 20069  8196
Mo_E18_5_MesCells <-  readRDS("Mo_E18_5_MesCells.rds")
dim(Mo_E18_5_MesCells)
## [1] 20069  7585

1.3.2 Normalizing the data

After removing unwanted cells from the dataset, the next step is to normalize the data. By default, we employ a global-scaling normalization method “LogNormalize” that normalizes the feature expression measurements for each cell by the total expression, multiplies this by a scale factor (10,000 by default), and log-transforms the result. Normalized values are stored in Mo_E18_5_MesCells[[“RNA”]]@data.

Mo_E18_5_MesCells <- NormalizeData(Mo_E18_5_MesCells, normalization.method = "LogNormalize",
                         scale.factor = 10000)

1.3.3 Identification of highly variable features (feature selection)

We next calculate a subset of features that exhibit high cell-to-cell variation in the dataset (i.e, they are highly expressed in some cells, and lowly expressed in others). We and others have found that focusing on these genes in downstream analysis helps to highlight biological signal in single-cell datasets.

Our procedure in Seurat3 is described in detail here, and improves on previous versions by directly modeling the mean-variance relationship inherent in single-cell data, and is implemented in the FindVariableFeatures function. By default, we return 2,000 features per dataset. These will be used in downstream analysis, like PCA.

Mo_E18_5_MesCells <- FindVariableFeatures(Mo_E18_5_MesCells, selection.method = "vst", nfeatures = 2000)
# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(Mo_E18_5_MesCells), 10)

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(Mo_E18_5_MesCells)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot2

1.3.4 Scaling the data

Next, we apply a linear transformation (‘scaling’) that is a standard pre-processing step prior to dimensional reduction techniques like PCA. The ScaleData function:

Shifts the expression of each gene, so that the mean expression across cells is 0 Scales the expression of each gene, so that the variance across cells is 1 This step gives equal weight in downstream analyses, so that highly-expressed genes do not dominate The results of this are stored in Mo_E18_5_MesCells[[“RNA”]]@scale.data

all.genes <- rownames(Mo_E18_5_MesCells)
Mo_E18_5_MesCells <- ScaleData(Mo_E18_5_MesCells, features = all.genes)

1.3.5 Perform linear dimensional reduction

Next we perform PCA on the scaled data. By default, only the previously determined variable features are used as input, but can be defined using features argument if you wish to choose a different subset.

Mo_E18_5_MesCells <- RunPCA(Mo_E18_5_MesCells, features = VariableFeatures(object = Mo_E18_5_MesCells))

Seurat provides several useful ways of visualizing both cells and features that define the PCA, including VizDimReduction, DimPlot, and DimHeatmap

# Examine and visualize PCA results a few different ways
print(Mo_E18_5_MesCells[["pca"]], dims = 1:5, nfeatures = 5)
## PC_ 1 
## Positive:  Serpinh1, Col6a1, Fstl1, Sparc, Col3a1 
## Negative:  Lyz2, Pglyrp1, Mmp8, Stfa3, Ifitm6 
## PC_ 2 
## Positive:  Stfa2l1, Cstdc4, Stfa2, Cstdc5, Csta3 
## Negative:  Pclaf, Top2a, Stmn1, Birc5, Cdca8 
## PC_ 3 
## Positive:  Ctsb, Ctss, Ms4a6c, Apoe, Ms4a6b 
## Negative:  Stfa2, Cstdc5, Retnlg, Csta3, Cd24a 
## PC_ 4 
## Positive:  Akap12, Entpd2, Igf1, Smoc2, Agtr2 
## Negative:  Hapln1, Sox5, Fam180a, Col9a1, Col2a1 
## PC_ 5 
## Positive:  Tspan15, Cd55, Cdh13, Ackr2, Htra4 
## Negative:  Col9a3, Col2a1, Col9a2, Col11a2, Col9a1
VizDimLoadings(Mo_E18_5_MesCells, dims = 1:2, reduction = "pca")

DimPlot(Mo_E18_5_MesCells,group.by = "orig.ident", reduction = "pca")

DimHeatmap(Mo_E18_5_MesCells, dims = 1:9, cells = 500, balanced = TRUE)

1.4 Determine the ‘dimensionality’ of the dataset

ElbowPlot(Mo_E18_5_MesCells, ndims = 50)

1.4.1 Cluster the cells

Seurat v3 applies a graph-based clustering approach, building upon initial strategies in (Macosko et al). Importantly, the distance metric which drives the clustering analysis (based on previously identified PCs) remains the same. However, the approach to partioning the cellular distance matrix into clusters has dramatically improved. The approach was heavily inspired by recent manuscripts which applied graph-based clustering approaches to scRNA-seq data [SNN-Cliq, Xu and Su, Bioinformatics, 2015] and CyTOF data [PhenoGraph, Levine et al., Cell, 2015]. Briefly, these methods embed cells in a graph structure - for example a K-nearest neighbor (KNN) graph, with edges drawn between cells with similar feature expression patterns, and then attempt to partition this graph into highly interconnected ‘quasi-cliques’ or ‘communities’.

As in PhenoGraph, we first construct a KNN graph based on the euclidean distance in PCA space, and refine the edge weights between any two cells based on the shared overlap in their local neighborhoods (Jaccard similarity). This step is performed using the FindNeighbors function, and takes as input the previously defined dimensionality of the dataset (first 25 PCs).

To cluster the cells, we next apply modularity optimization techniques such as the Louvain algorithm (default) or SLM [SLM, Blondel et al., Journal of Statistical Mechanics], to iteratively group cells together, with the goal of optimizing the standard modularity function. The FindClusters function implements this procedure, and contains a resolution parameter that sets the ‘granularity’ of the downstream clustering, with increased values leading to a greater number of clusters. We find that setting this parameter between 0.4-1.2 typically returns good results for single-cell datasets of around 3K cells. Optimal resolution often increases for larger datasets. The clusters can be found using the Idents function.

Mo_E18_5_MesCells <- FindNeighbors(Mo_E18_5_MesCells, dims = 1:25)
Mo_E18_5_MesCells <- FindClusters(Mo_E18_5_MesCells, resolution = 0.1)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 7585
## Number of edges: 283352
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9755
## Number of communities: 11
## Elapsed time: 0 seconds
Mo_E18_5_MesCells <- FindClusters(Mo_E18_5_MesCells, resolution = 0.5)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 7585
## Number of edges: 283352
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9277
## Number of communities: 22
## Elapsed time: 0 seconds
Mo_E18_5_MesCells <- FindClusters(Mo_E18_5_MesCells, resolution = 0.2)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 7585
## Number of edges: 283352
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9594
## Number of communities: 15
## Elapsed time: 0 seconds
Mo_E18_5_MesCells <- FindClusters(Mo_E18_5_MesCells, resolution = 0.3)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 7585
## Number of edges: 283352
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9484
## Number of communities: 18
## Elapsed time: 0 seconds

1.4.2 Run non-linear dimensional reduction (UMAP/tSNE)

Seurat offers several non-linear dimensional reduction techniques, such as tSNE and UMAP, to visualize and explore these datasets. The goal of these algorithms is to learn the underlying manifold of the data in order to place similar cells together in low-dimensional space. Cells within the graph-based clusters determined above should co-localize on these dimension reduction plots. As input to the UMAP and tSNE, we suggest using the same PCs as input to the clustering analysis.

# If you haven't installed UMAP, you can do so via reticulate::py_install(packages =
# 'umap-learn')
Mo_E18_5_MesCells <- RunUMAP(Mo_E18_5_MesCells, dims = 1:25)
# note that you can set `label = TRUE` or use the LabelClusters function to help label
# individual clusters
DimPlot(Mo_E18_5_MesCells, reduction = "umap", label = T, label.size = 5)

# note that you can set `label = TRUE` or use the LabelClusters function to help label
# individual clusters
FeaturePlot(Mo_E18_5_MesCells, reduction = "umap",
            features = c("nCount_RNA", "nFeature_RNA"))

1.5 Find Doublets

I use here DoubletFinder for predcting doublets in the dataset

DoubletFinder can be broken up into 4 steps: (1) Generate artificial doublets from existing scRNA-seq data (2) Pre-process merged real-artificial data (3) Perform PCA and use the PC distance matrix to find each cell’s proportion of artificial k nearest neighbors (pANN) (4) Rank order and threshold pANN values according to the expected number of doublets

## pK Identification (no ground-truth) ---------------------------------------------------------------------------------------
sweep.res.list_Fpl <- paramSweep_v3(Mo_E18_5_MesCells, PCs = 1:10, sct = FALSE)
sweep.stats_Fpl <- summarizeSweep(sweep.res.list_Fpl, GT = FALSE)
bcmvn_Fpl <- find.pK(sweep.stats_Fpl)
## Homotypic Doublet Proportion Estimate -------------------------------------------------------------------------------------
homotypic.prop <- modelHomotypic(Mo_E18_5_MesCells@meta.data$RNA_snn_res.0.2)           ## ex: annotations <- seu_kidney@meta.data$ClusteringResults
nExp_poi <- round(0.075*nrow(Mo_E18_5_MesCells@meta.data))  ## Assuming 7.5% doublet formation rate - tailor for your dataset
nExp_poi.adj <- round(nExp_poi*(1-homotypic.prop))
Mo_E18_5_MesCells <- doubletFinder_v3(Mo_E18_5_MesCells, PCs = 1:10,
                            pN = 0.25, pK = 0.09,
                            nExp = nExp_poi, reuse.pANN = FALSE,
                            sct = FALSE)
Mo_E18_5_MesCells <- doubletFinder_v3(Mo_E18_5_MesCells, PCs = 1:10, 
                            pN = 0.25, pK = 0.09,
                            nExp = nExp_poi.adj,
                            reuse.pANN = "pANN_0.25_0.09_569",
                            sct = FALSE)
DimPlot(Mo_E18_5_MesCells, group.by="DF.classifications_0.25_0.09_569",
        cols = c("red","grey"),
        reduction="umap", pt.size=0.3 
        )

1.6 Save Object for further analyses

saveRDS(Mo_E18_5_MesCells, "Mo_E18_5_MesCells.rds")

2 Markers

From Gretel:

There are some “pan” markers for the different compartments and clusters and highly expressed markers (altough not exclusive) for clusters already well defined in the literature, I send you some examples and explanations.

These clusters were defined following the unsupervised clustering first and redifined in more subclusters from our knowledge of the literature. Without doublets or dead cell filtration.

In the stromal compartment there are clusters differents from adult and we have many questions, the aim of the study.

2.1 Hematopoietic compartment

2.1.1 a) Lymphoid and myeloid defined by CD45/ptprc

Other markers are: CD52, CD53

FeaturePlot(Mo_E18_5_MesCells,
            features = c("Ptprc", # Leukocyte lineage
                         "Cd52", # 
                         "Cd53") # 
            )

2.1.2 B) Erythroid lineage

Is identified by Hbb-bt, Hba-a1, Hba-a2 (by FACS with TER119)

FeaturePlot(Mo_E18_5_MesCells,
            features = c("Hbb-bt", # 
                         "Hba-a1", # 
                         "Hba-a2") # 
            )

2.2 2) Endothelial compartment:

Tek, Tie1, Cdh5, Sox7, Pecam1, (Eng, CD200, Ly6a, not exclusives but highly expressed in these cells)

In perinatal is only one cluster but in adult 2 clusters are well defined: the arteriolar and the sinusoidal endothelial cells.

FeaturePlot(Mo_E18_5_MesCells,
            features = c("Tek", # 
                         "Tie1", # 
                         "Cdh5",
                         "Pecam1") # 
            )

2.3 3) Mesenchymal Compartment

Pan markers: Pdgfra, Prrx1, Col1a2,

FeaturePlot(Mo_E18_5_MesCells,
            features = c("Pdgfra", # 
                         "Prrx1", # 
                         "Col1a2") # 
            )

Committed lineages are easier to identify.

2.3.1 Osteoblasts & Hyperthrophic chondrocytes (HC is mainly in perinatal stages):

Sp7, Runx2, Bglap, Bglap2,Bglap3, Mmp13, Col13a, Ibsp, Col1a1 (in perinatal stages they are all together) but are mainly in two populations in adult (according to the literature are progenitors of the osteo lineage).

FeaturePlot(Mo_E18_5_MesCells,
            features = c("Sp7", # 
                         "Mmp13", # 
                         "Bglap",
                         "Bglap2")# 
            )

2.3.2 Chondroblasts and chondrocytes:

Sox9,Sox5,Sox6, Col2a1, Acan, Col9a3, Col9a1,Nt5e (In perinatal stages some of the markers are not exclusive but the combination of all define the cluster).

FeaturePlot(Mo_E18_5_MesCells,order = T,
            features = c("Sox9", # 
                         "Sox5", # 
                         "Sox6",
                         "Col2a1")# 
            )

Myofibroblasts: Pax7, Myf5, Myod1, Mrln, Msc

FeaturePlot(Mo_E18_5_MesCells,order = T,
            features = c("Pax7", # 
                         "Myf5", # 
                         "Myod1",
                         "Mrln")# 
            )

Other clusters in the bone

Neural cells: Schwann cells: Sox10, Mal, Plp1

FeaturePlot(Mo_E18_5_MesCells,order = T,
            features = c("Sox10", # 
                         "Mal", # 
                         "Plp1" )
            )

DimPlot(Mo_E18_5_MesCells, label = T, group.by = "RNA_snn_res.0.3" )

dittoBarPlot(object = Mo_E18_5_MesCells,
             var = "DF.classifications_0.25_0.09_569",
             group.by = "RNA_snn_res.0.3")

2.4 Export to rds for viz

source(file = "../utils/seurat2shiny.R")

seurat2shiny(object = Mo_E18_5_MesCells, assay = "RNA",save = T,
             name = "rds_viz/Mo_E18_5_MesCells_preprocess")

2.5 Make a list of gene markers for each cluster

all.markers <- FindAllMarkers(object = Mo_E18_5_MesCells, only.pos = T, 
                              logfc.threshold = 0.5, min.diff.pct = 0.2)

write.csv(x = all.markers, file ="Mo_E18_5_MesCells_allmarkers.csv")
sessionInfo()
## R version 4.1.0 (2021-05-18)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=es_ES.UTF-8        LC_COLLATE=en_US.UTF-8     LC_MONETARY=es_ES.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=es_ES.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=es_ES.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] RColorBrewer_1.1-2  dittoSeq_1.5.1      DoubletFinder_2.0.3 ggplot2_3.3.5       patchwork_1.1.1     SeuratObject_4.0.2  Seurat_4.0.3        dplyr_1.0.7         Matrix_1.3-4        BiocStyle_2.20.2   
## 
## loaded via a namespace (and not attached):
##   [1] plyr_1.8.6                  igraph_1.2.6                lazyeval_0.2.2              splines_4.1.0               listenv_0.8.0               scattermore_0.7             GenomeInfoDb_1.28.1         digest_0.6.27               htmltools_0.5.1.1           magick_2.7.2                fansi_0.5.0                 magrittr_2.0.1              tensor_1.5                  cluster_2.1.2               ROCR_1.0-11                 globals_0.14.0              matrixStats_0.59.0          spatstat.sparse_2.0-0       colorspace_2.0-2            ggrepel_0.9.1               xfun_0.24                   crayon_1.4.1                RCurl_1.98-1.3              jsonlite_1.7.2              spatstat.data_2.1-0         survival_3.2-11             zoo_1.8-9                   glue_1.4.2                  polyclip_1.10-0             gtable_0.3.0                zlibbioc_1.38.0            
##  [32] XVector_0.32.0              leiden_0.3.8                DelayedArray_0.18.0         future.apply_1.7.0          SingleCellExperiment_1.14.1 BiocGenerics_0.38.0         abind_1.4-5                 scales_1.1.1                pheatmap_1.0.12             DBI_1.1.1                   miniUI_0.1.1.1              Rcpp_1.0.7                  viridisLite_0.4.0           xtable_1.8-4                reticulate_1.20             spatstat.core_2.2-0         stats4_4.1.0                htmlwidgets_1.5.3           httr_1.4.2                  ellipsis_0.3.2              ica_1.0-2                   pkgconfig_2.0.3             farver_2.1.0                sass_0.4.0                  uwot_0.1.10                 deldir_0.2-10               utf8_1.2.1                  tidyselect_1.1.1            labeling_0.4.2              rlang_0.4.11                reshape2_1.4.4             
##  [63] later_1.2.0                 munsell_0.5.0               tools_4.1.0                 generics_0.1.0              ggridges_0.5.3              evaluate_0.14               stringr_1.4.0               fastmap_1.1.0               yaml_2.2.1                  goftest_1.2-2               knitr_1.33                  fitdistrplus_1.1-5          purrr_0.3.4                 RANN_2.6.1                  pbapply_1.4-3               future_1.21.0               nlme_3.1-152                mime_0.11                   compiler_4.1.0              plotly_4.9.4.1              png_0.1-7                   spatstat.utils_2.2-0        tibble_3.1.2                bslib_0.2.5.1               stringi_1.6.2               highr_0.9                   lattice_0.20-44             vctrs_0.3.8                 pillar_1.6.1                lifecycle_1.0.0             BiocManager_1.30.16        
##  [94] spatstat.geom_2.2-0         lmtest_0.9-38               jquerylib_0.1.4             RcppAnnoy_0.0.18            data.table_1.14.0           cowplot_1.1.1               bitops_1.0-7                irlba_2.3.3                 httpuv_1.6.1                GenomicRanges_1.44.0        R6_2.5.0                    bookdown_0.22               promises_1.2.0.1            KernSmooth_2.23-20          gridExtra_2.3               IRanges_2.26.0              parallelly_1.26.1           codetools_0.2-18            MASS_7.3-54                 assertthat_0.2.1            SummarizedExperiment_1.22.0 withr_2.4.2                 sctransform_0.3.2           S4Vectors_0.30.0            GenomeInfoDbData_1.2.6      mgcv_1.8-36                 parallel_4.1.0              grid_4.1.0                  rpart_4.1-15                tidyr_1.1.3                 rmarkdown_2.9              
## [125] MatrixGenerics_1.4.0        Rtsne_0.15                  Biobase_2.52.0              shiny_1.6.0